ON 
On 

>v A two-dimensional network simulator for two-phase 

& : 

flow in porous media 

m 

CN ■ Eyvind Aker* Knut J0rgen Mal0y 

Department of Physics 
Ch ; University of Oslo, N-0316 Oslo, Norway 

: 

Alex Hansen* 

y- j — — | Department of Physics, Norwegian University of Science 

^ and Technology, N-7034 Trondheim, Norway 

O ' 

c/2 . G. George Batrouni 

Institut Non-Lineaire de Nice, 
Q H ' Universite de Nice - Sophia Antipolis, 06560 Valbonne, France 

February 2, 2008 

> 

m 
o 

I/") . Abstract 

o 

' We investigate a two-dimensional network simulator capable of 

0^ . modeling different time dependencies in two-phase drainage displace- 

""^5 ' ments. In particular, we focus on the temporal evolution of the pres- 

O , sure due to capillary and viscous forces and the time dependence of 

c/) ■ the interface between the two liquids. The dynamics of the capillary 

effect are taken into account and we report on high accuracy pressure 
measurements. Moreover, the simulator includes important features in 
drainage, like burst dynamics of the invading fluid and simultaneous 
flow of two liquids in a section of a tube. The validity of the model is 
checked by comparing simulation results with well known experimental 
'jj] ■ properties in drainage displacement. 
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1 Introduction 



Two-phase displacements in porous media have been much studied over the 
last two decades. The main reason for this is the great variety of structures 
observed when changing the physical parameters of the fluids like viscos- 
ity contrast, wettability, interfacial tension and displacement rate. Besides 
being a process of great interest in modern physics, it has a large num- 
ber of practical applications in many fields of science like oil recovery and 
hydrology. 

The main purpose of this paper is to present a network simulator mod- 
eling immiscible two-phase flow on a two-dimensional lattice of cylindrical 
tubes. Primarily, the model is developed to measure the time dependence 
of different physical properties and to study the dynamics of the fluid move- 
ments. We focus on drainage displacements, i.e. the process where a non- 
wetting fluid displaces a wetting fluid in a porous medium. In particular, we 
report on the dynamics of the temporal evolution of the pressure due to cap- 
illary and viscous forces and the time dependence of the interface between 
the two liquids. 

The different structures between the invading and the defending fluids 
obtained in drainage divide into three major flow regimes: viscous finger- 
ing H,16|, stable displacement 13] and capillary fingering [15|. There exist 
statistical models like DLA |22( |, anti-DLA |18| and invasion percolation [ j21| j 
that reproduce the basic domains in viscous fingering, stable displacement 
and capillary fingering respectively. However, these models do not contain 
any physical time for the front evolution and the they cannot describe the 
crossover between the major flow regimes. 

To overcome the limitations of the statistical models several network 
simulators similar to the one presented in this paper, have been developed 
over the last two decades P,^-|6|,pT|-|T3|. The different models are more or less 
simplified to avoid computational complications, and difficulties most often 
arise when the capillary effects are taken into account. Lenormand et al. 
(1988) assigned each tube a threshold pressure and allowed the injecting non- 
wetting fluid to enter the tube only when the pressure exceeded the threshold 
value of that tube. Moreover, the non-wetting fluid was restricted to flow in 
the positive direction relative to the direction of the displacement. Dias and 
Payatakes (1986) suggested a more realistic approximation by introducing 
tubes with walls with a sinusoidal profile. They let the capillary pressure of 
the meniscus change as the menisci invaded the tubes. 

Similar to the idea of Dias and Payatakes (1986) the model reported here 
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takes into account the capillary effects by letting the capillary pressure of 
the meniscus depend on its position inside the tube. Moreover, the menisci 
can both invade into or retreat from a tube. This is an important property 
in slow drainage, where the invading fluid is found to suddenly invade a 
larger region causing another part to retire (bursts) 1(| 14], 17|. Experiments 
performed by Mal0y et al. (1992) have also produced evidence that the bursts 
are characterized by large pressure fluctuations (Haines jumps). To measure 
these pressure fluctuations the flow field in our network simulator is solved 
for a constant injection rate rather than a constant pressure. 

In the present model an approximation is developed to model a mixture 
of the wetting and non-wetting fluids when they flow simultaneously in a 
cross section of a tube. Such mixing is an important process in both drainage 
and imbibition f?,14] and should not be neglected. The piston-like motion 
characterizing slow drainage is well described, but additional mechanisms 
observed in imbibition |]3, 14 1 is not taken into account. 

The paper is organized as follows. Section 2 describes the porous medium 
model used in the network simulator. Section 3 gives the solution of the 
flow field with the constraint of a constant injection rate, and Section 4 
describes the algorithm for updating the menisci. In Section 5 we discuss 
how to move the menisci into neighboring tubes and allow mixing of the 
two liquids. Finally, in Section 6 we give some simulation result focusing on 
calibrations of the model. At the end some conclusions are discussed. 



2 The Porous Medium Model 

2.1 Geometry of Model Porous Medium 

The porous medium is represented by a square lattice of tubes inclined 45 
degrees. Thus, if all tubes are equal and a uniform pressure across the 
lattice is applied, a liquid flows equally well in tubes inclined to the left 
as tubes inclined to the right. The tubes are connected together at nodes, 
where four tubes meet. There is no volume assigned to the nodes: the 
tubes represent the volume of both pores and throats. The liquids flow 
from the bottom to the top of the lattice and periodic boundary conditions 
are applied horizontally. The pressure difference between the first (bottom 
of system) and the last (top) rows defines the pressure across the lattice. 
Gravity effects are neglected, and as a consequence we consider a horizontal 
flow in a two-dimensional network of tubes. See figure |l| for details. 

The tubes are cylindrical with equal length d, and every tube is assigned a 
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radius r which is chosen at random in the interval [Ai, \?\d. The randomness 
of the radii represent the disorder of an ordinary porous medium and Ai and 
A2 define the width of the distribution of radii. 

Initially, the system is filled with a defending fluid with viscosity n\. 
The invading fluid with viscosity \i2 is injected from the bottom row with a 
constant injection rate. Here we report on the study of drainage displace- 
ments, and so let the invading fluid be non-wetting and the defending fluid 
be wetting. We assume that the fluids are immiscible and that there is a well 
defined interface between the two phases. The curvature of this interface 
gives rise to a capillary pressure given by the interfacial tension 7. Moreover, 
we will treat the liquids as incompressible, i.e. the volume flux flowing into 
the system must equal the volume flux flowing out of the system. 



2.2 Capillary Flow in a tube 

Consider a tube containing a meniscus between nodes i and j in the network 
as shown in figure ||. The volume flux q^ through the tube from the ith to 
the jth node is found from the Washburn equation for capillary flow [pC|]: 

Vij = — " ■ l( A Pij ~ Pc) • (!) 

Mf d 

Here kij = rfj/8 is the permeability, fi e ff = fj,2Xij+fii(l — Xij) is the effective 
viscosity due to the two fluids, Apij = pj — is the pressure difference 
between the ith. and jth node and p c = p\—p2 defines the capillary pressure. 
Xij is the position of the meniscus in the tube and rij is the radius of the 
tube. The position of the meniscus is a continuous function in the range [0, 1] 
and the flow direction is given by the sign of qif When q^ > the fluids in 
the tube flow to the right, otherwise they flow to the left. As mentioned in 
the introduction, it is important to allow flow in both directions in order to 
study the dynamics when the menisci invade into or retreat from a tube. 

For a tube without a meniscus present p c = 0, and equation ([j]) reduces 
to that describing Hagen-Poiseulle flow. 

The capillary pressure p c due to the meniscus is given by the Young- 
Laplace law. Let 6 refer to the wetting angle between the non-wetting and 
wetting phases as shown in figure ||. Then 

27 

Pc = — cos 9 , (2) 
r 
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where rj cos 9 is the principal radius of curvature of the meniscus. 7 denotes 
the interfacial tension between the two phases. 

Equation (|2|) is derived under the assumption that the fluids are in static 
equilibrium, i.e. Apij = p c . A more detailed description of the dynamics of 
a moving interface is still an open problem (8|. However, at low flow rates 
when no mixing or turbulence occur, we consider the above expression as a 
reasonable approximation. 

A real porous medium consists of a complex network of throats and 
pores with a great variety of shapes. The curvature of a meniscus traveling 
in this network, becomes a continuous function depending on the meniscus' 
position. The variation in curvature results in local changes in the capillary 
pressure. Equation @ does not take this effect into account, since the 
capillary pressure is assumed constant for each tube. Instead, we apply a 
dependency in p c as a function of the menisci's position inside each tube. 
Thus, we define the capillary pressure p c as (without index notations): 

Pc = [1 - cos(2vrx)] , (3) 

where we have assumed perfect wetting (6 = 0). In the above formula r 
denotes the radius of the actual tube and x is the position of the meniscus 
in that tube. The function is plotted in figure |3[ The definition sets the 
capillary pressure equal to zero at the ends of the tube whereas p c becomes 
equal to the threshold pressure pt when the meniscus is in the middle of 
the tube, i.e. pt = 4rf/r. The threshold pressure is the minimum capillary 
pressure required to let the non-wetting fluid invade the tube. 

3 Solving the Flow Field 

The fluids are assumed incompressible leading to conservation of volume 
flux at each node: 

E% = °- ( 4 ) 

3 

Here qij denotes the flow trough a tube connecting node i and j. Equa- 
tion (Ej) is simply Kirchhoff's equation and in the summation j runs over 
the nearest neighbor nodes to the ith node. The index i runs over all nodes 
that do not belong to the top or bottom rows, that is, the internal nodes. 
This set of linear equations are to be solved with the constraint that the 
pressures at the nodes belonging to the upper and lower rows are kept fixed. 
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By inserting equation (|l]) in equation (Q) we get: 

E Gi d & ~ Pi ~ Pc) = , (5) 

3 

where Gy defines the mobility of the tube, i.e. Gij = Trrfjkij/fJ, e ff. In 
order to write this set of equations as a matrix equation, we move all the 
capillary pressures and the fixed pressure referring to the nodes belonging 
to the upper and lower rows to the right-hand side of equation (||). The 
final matrix equation may then be written as 

Yl Di iPi = Bi ' ( 6 ) 

3 

where the indices i and j only run over internal nodes. are elements 
a conductance matrix D [jl] where the elements depend on the connection 
between different tubes and their respective mobility, pj are the elements 
in the pressure vector, containing the pressure at the internal nodes and 
Bi contains the pressure at the boundaries (upper and lower rows) and the 
capillary pressure if a meniscus is present in the tube. 

We are seeking the solution of the pressure at the internal nodes for a 
given configuration of the menisci, i.e 

P 3 = E i D '%Bi • (7) 

i 

This equation is solved by using the Conjugate Gradient method 
3.1 Solving for a Constant Injection Rate 

The pressures solved by equation (0) correspond to keeping the applied pres- 
sure difference across the network constant. We want to study the dynamics 
of the pressure fluctuations at a constant displacement rate. Thus, we need 
to find the pressures pj under constant injection rate of the invading fluid. 

For two-phase displacement in a porous medium the injection rate Q is 
given by 

Q = AAP + B . (8) 

Here AP is the pressure across the lattice and A and B are parameters 
depending on the geometry of the medium and the current configuration of 
the liquids. The first part of equation (|8|) is simply Darcy's law for one phase 
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flow through a porous medium. The last part B results from the capillary 
pressure between the two phases. As long as the menisci are kept at fixed 
positions B becomes a constant. 

The pressure across the lattice for a given injection rate is 

AP=^. (9) 

The parameters A and B are calculated by solving the pressure field (^) 
for two different pressure AP and AP applied across the lattice. The 
obtained pressure vectors give the corresponding injection rates, Q and Q . 
We are now in the position to calculate A and B by solving the set 

Q' = AAP' + B , (10) 
Q" = AAP" + B . (11) 

The next step is to relate the internal pressures at the nodes to the 
constant injection rate. From equation (|l|) we can write the flow rate q^ 
from node i to node j like 

qij = aijApij + bij , (12) 

where a,ij and &y are parameters depending on the permeability of the tube, 
the effective viscosity and the capillary pressure due to a meniscus. From 
equation (^) the pressure AP for a given Q is found. Thus, we could proceed 
by solving Kirchhoff 's equation with the constraint of applying this pressure 
at the inlet. That would give the correct flow rates in each tube for the 
desired Q. However, the following method will save a third solution of 
equation (0). 

All the equations involved in the calculations are linear, i.e. they have 
the functional form f(x) = ax + b. As a consequence the pressure Apij 
between node i and j becomes a linear function of the pressure AP across 
the lattice: 

Ap^ = FijAP + Uij , (13) 

Here are Tjj and n^- parameters depending on the configuration of the fluids. 
By inserting this into equation ( |l2"|) and redefine ctf/'s and fry's we obtain 

qij = aijAP + kj . (14) 

The parameters Sy and bij are found from the flow rates q[j and q'/j cor- 
responding to the two pressures AP' and AP" respectively. Note that the 
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parameters A and B in (||) and S^- and fry in ( |14| ) all depend on the current 
position of the menisci and we therefore need to solve them for every new 
fluid configuration. 

The solution due to a constant injection rate can now be summarized 
into two steps: (a) After we have found A and B we use equation (^) to get 
AP for the desired Q. (b) This AP is then used in equation (|l4|) to get the 
local flow qij. 

The validity of equation ( |T3| ) is easily checked by solving Kirchhoff's 
equation for the calculated pressure AP and compare the solution with the 
one given from equation (|14|). Numerical results show excellent agreement 
between these two solutions. 



4 Updating the Menisci 

A time step At is chosen such that every meniscus is allowed to travel at 
most a maximum step length Ax max during that time step. This leads to 
the formula 

Ax inn 



At = min 

ij 



(15) 



where = qij/^rf- denotes the flow velocity in a tube containing a meniscus 
between the ith and the jth node. The time step becomes dependent on the 
local velocity. This method is sometimes called event driven updating. 

During the time step it is checked whether or not a meniscus crosses 
the middle or the end of a tube. If this happens, the time step is redefined 
such that only one meniscus reaches the end or the middle of the tube. A 
meniscus reaching the end of a tube is moved into the neighbor tubes (see 
section |||) . 

In the middle of the tube the capillary pressure becomes equal to the 
threshold pressure. At this point, the menisci in slow drainage are found 
to become unstable and suddenly invade the tube like a burst 0, [L4|, [lTj . 
Often a cascade of bursts are released in rapid succession around the original 
instability, before a new stable position is reached. Moreover, the bursts are 
characterized by large pressure fluctuations [17|. It is important to detect 



the order in which the instabilities occur to ensure that the right path of 
least resistance is chosen. This is done by adjusting the time step when a 
meniscus approaches the middle of a tube such that the exact occurrence of 
the burst is detected. 
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The new positions of the menisci are calculated by using a second order 
Runge-Kutta scheme |jT9f ■ Numerical analysis from the present model shows 
that this scheme produces more stable solutions than the less accurate and 
more unstable Euler scheme. However, with an event driven time step as 
defined in equation (|l5|), the implementation of the second order Runge- 
Kutta scheme is not straight forward. Before discussing our modifications, 
we first remind the reader of the Runge-Kutta scheme. 

Let x n denote the position of a meniscus at time t n and let the next 
position at time t n+ \ = t n + At be x n+ \ where At is the time step found 
from (15). The second order Runge-Kutta scheme then becomes [19| 



x n+ i = x n + k 2 + 0(At 3 ) , (16) 
k\ = At ■ v(t n ,x n ) , (17) 
k 2 = At-v(t n + \At,x n + \ki) , (18) 

where v(t n , x n ) denotes the local flow rate in each tube and v(t n + ^ At, x n + 
— Vmiditni %ri) defines the midpoint velocity. 
Note that the velocities defining the time step in equation (|l~5| ) corre- 
spond to the derivative of the curve x n at the starting point of each time 
interval. I.e. v(t n ,x n ) = Vij, where the subscript ij is omitted on the left 
hand side of the equality. When the menisci are updated by using the ve- 
locities v(t n ,x n ), equation (FL5l) leads to the condition 



l^n+l %n\ 5; Ax maa; , (19) 

for all the menisci. In the second order Runge-Kutta scheme the next po- 
sition x n+ \ is found by using the midpoint velocity v rn q i( i(t n ,x n )^ which in 
general differs from v(t n ,x n ). As a consequence, condition ( |l~9| ) may no 
longer be valid and the displacement length or the stability of the numerical 
solution, are no longer controlled. Assume that the position of the menisci is 
a smooth function of time, the effect vanishes since v m id(t n ,x n ) ~ v(t n ,x n ). 
Moreover, there is no problem as long as v m idkPni X n ) < v(t n ,x n ). The prob- 
lem arises when v m id{t n ,x n ) ^> v(t n ,x n ) and the final displacement becomes 
much larger than Ax max . To avoid this scenario the time step is redefined 
by inserting the midpoint velocities in equation (|l~5|). A new midpoint veloc- 
ity corresponding to the redefined time step is calculated and the positions 
are updated according to this time step by using the new midpoint velocity. 
The procedure is repeated until the displacements become close enough to 
the maximum step length Ax max . Typically, the final displacement must 
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not be larger than about 10% increase of Ax max . Otherwise, the numerical 
solution may diverge. 

5 Motion of the Menisci at the Nodes 

A tube partially filled with one of the liquids is allowed to have either one or 
two menisci leading to four different arrangements as shown in figure ||. The 
menisci can be situated at any position inside the tube and the correspond- 
ing effective viscosity and capillary pressure are calculated. The effective 
viscosity becomes the sum of the fraction of the wetting and non-wetting 
fluid inside the tube multiplied with their respective viscosities. The abso- 
lute value of the capillary pressure is given by equation (||), while its sign 
depends on whether the meniscus is pointing upwards like in figure || (a) or 
downwards like in figure |(b). The menisci are updated according to the 
determined time step At from the previous section, and their respective flow 
rates Vij. The total time lapse is recorded before the flow field is solved for 
the new fluid configuration. Note that due to the incompressibility of the 
liquids the two menisci within the same tube in figure |] (c) and (d) always 
move with equal velocities. Thus, the volumes of the wetting in (c) or the 
non-wetting in (d) are conserved. 

When a meniscus reaches the end of a tube it is moved into the neighbor 
tubes according to some defined rules. These rules take care of the different 
fluid arrangements that can appear around the node. Basically, the non- 
wetting fluid can either invade into or retreat from the neighbor tubes as 
shown in figure || (a) and (b) respectively. In ||| (a) the non- wetting fluid 
approaches the node from below (drainage). When the meniscus has reached 
the end of the tube (position 1), it is removed and three new menisci are 
created at position 5 in the neighbor tubes (position 2). The distance 5 is 
about 1-5% of the tube length d and it defines the node region where the 
capillary pressure is zero. The node region avoids that the created menisci 
at position 2 immediately disappear and move back to the initial position 1 
in tubes where the flow direction is opposite to the direction of the invading 
fluid. The node region has also an important role that allows mixing of the 
liquids (see below). 

Figure [| (b) shows the opposite case when the non-wetting fluid retreats 
into a single tube (imbibition) . As figure || shows the properties of imbibition 
should not be neglected as long as the menisci can travel in both directions. 
However, in drainage which is what we are focusing on, arrangement (b) will 
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appear rarely compared to (a). 

When the menisci are moved a distance 5 into the neighbor tubes the 
total time lapse is adjusted due to the injection rate of the invading fluid. 
The adjustment is calibrated such that the amount of the invading fluid in 
the lattice always equals the injected volume. The injected volume is the 
product of the time lapse and the injection rate. 

Difficulties arise when we want to create a new meniscus in a tube that 
already contains two menisci. To allow such movement the two original 
menisci and the new meniscus are merged into one by the following method. 
Consider a scenario like figure [| (a) but assume now that the leftmost tube 
already contains a bubble of non- wetting fluid as shown in figure ^(a). To 
merge the three menisci in the leftmost tube into one, we move the wetting 
fluid between position 2 and 3 to the left side of the non-wetting bubble and 
remove the menisci at position 2 and 3 before the meniscus at position 4 is 
moved to the right a distance equal to the original length between position 
2 and 3 (figure |6](b)). The same principles apply when the non- wetting fluid 
retreat into a tube that already contains a bubble of wetting fluid. 

5.1 Mixing of the Fluids 

The moving rules described in figure [| and || solve the problem of modeling 
a "mixture" of the two liquids. Consider the situations shown in figure |7] 
where both the non-wetting and wetting fluids flow toward the node from the 
bottom and right tube respectively. Physically, it is expected that the fluid 
flowing out of the node is a mixture of both liquids. It is observed fA [141] 
that such simultaneous flow of both liquids in a tube takes place in so- 
called funicular continuous flow or a discontinuous dispersed flow. In the 
former case the wetting phase flows along the cylinder wall surrounding 
the non-wetting phase, which occupies the central portion of the tube. A 
discontinuous flow is characterized by a dispersed non-wetting phase flowing 
as isolated droplets in a continuous wetting phase. 

A model describing funicular and dispersed flow would be too compli- 
cated. Instead, we assume that the simultaneous flow can be represented by 
a finite number of small bubbles of each liquid, placed next to each other 
inside the tube. By sorting the bubbles of same type of fluid they can be 
replaced by one or two menisci. The procedure is illustrated with an ex- 
ample in figure [?]: When the meniscus in the bottom tube reaches the end 
of the tube it moves a distance 5 into the neighbor tubes (arrangements 
a and b). Due to the opposite flow direction in the right tube the created 
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meniscus in this tube flows back to the node (arrangement c) and moves into 
the neighbor tubes (arrangement d). Now, the new meniscus in the bottom 
tube again approaches the node from below (arrangement e) and creates a 
configuration with three menisci in the top and left tubes (arrangement f). 
To avoid that the number of menisci inside a single tube increases unchecked 
the three menisci at the top and the left tube are reduced to one by placing 
the wetting fluid on the top of the non- wetting one (arrangement g). 

The reorganization when three menisci are reduced to one (figure ^ and ^) 
results in unphysical jumps in the capillary pressure. Due to the small size of 
the bubbles, the jumps usually appear as perturbations in the total pressure. 
A small number of such jumps will not affect the numerical solution very 
much. But, the moment they become more dominant the numerical solution 
may become unstable and diverge. We discuss this below. 



6 Calibration of the Network Model 

We present three different calibration simulations, one in each of the regimes 
of interest: viscous fingering, stable displacement and capillary fingering. 
The simulations are executed on a Cray T90 vector machine. The Con- 
jugate Gradient method solving Kirchhoff's equations, is easy to vectorize 
and achieves high performance efficiency on vector machines. However, the 
amount of time required increases dramatically with the size of the lattice. 
For a lattice consisting of N tubes, the Conjugate Gradient method needs 
approximately iV 2 iterations each time step to solve Kirchhoff's equations. 
Doubling the size of the lattice quadruples N, hence, the computation time 
increases by a factor of 16. In addition the number of the time steps before 
the invading fluid reaches the outlet strongly affects the CPU-time. 

In two-phase fluid displacement there are mainly three types of forces: 
viscous forces in the invading fluid, viscous forces in the defending fluid 
and capillary forces due to the interface between them. This leads to two 
dimensionless numbers that can characterize the flow in porous media: the 
capillary number C a and the viscosity ratio M. 

The capillary number is a quantity describing the competition between 
capillary and viscous forces. It is defined as 

Ca = ^ , (20) 

where Q (cm 2 /s) denotes the injection rate, /i (Poise) is the maximum vis- 
cosity of the two fluids , £ (cm 2 ) is the cross section of the inlet and 7 
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(dyn/cm) is the interfacial tension between the two phases. S is the prod- 
uct of the length of the inlet and the mean thickness of the lattice due to 
the average radius of the tubes. 

M defines the ratio of the viscosities of the two fluids and is given by 
the invading viscosity H2 divided with the defending viscosity fii- 

M = ^. 

Mi 

The three simulations are performed with parameters as close as possible 
to the corresponding experiments done by Mal0y et al. (1985) and Frette et 
al. (1997). In light of that, the length d of all tubes in the lattices are set 
equal to 1 mm. Furthermore, the radii r of the tubes are chosen randomly in 
the interval 0.05d < r < d. The interfacial tension is set to 7 = 30 dyn/cm 
and the viscosities of the defending and the invading fluids varies between 
0.01 P (~ water) and 10 P (~ glycerol). 

6.1 Viscous Fingering 

Figure |8| shows the result of a simulation in the regime of viscous finger- 
ing performed on a lattice of 60 x 80 nodes. The corresponding pressure 
across the lattice as a function of time is shown in figure ||| The simu- 
lation is stopped at breakthrough of the invading non-wetting fluid. The 
displacements are done with a high injection rate, Q = 1.5ml/min and the 
invading fluid is less viscous than the defending wetting fluid, ^2 = 0.010 P 
and hi = 10 P. The capillary number and the viscosity ratio becomes 
C a = 4.6 • 10~ 3 and M = 1.0 • 10~ 3 respectively. 

In viscous fingering the principal force is due to the viscous forces in 
the defending fluid and the capillary forces at the menisci are less dominant. 
The pattern formation (figure |8|) shows that the invading fluid creates typical 
fingers into the defending fluid. 

The pressure across the lattice (figure [|) decreases as the less viscous 
fluid invades the system. Roughly, the pressure appears to decrease linearly 
as a function of time. However, the slope is non-trivial and results from the 
fractal development of the fingers. In addition to the fractal growth, the 
rate of change in the pressure depends on the viscosity contrast M between 
the two phases. The rapid decrease at the end of the pressure function 
corresponds to the breakthrough of the invading fluid at the outlet. 

Figure |9] shows that there are small fluctuations in the average decreas- 
ing pressure function. The fluctuations correspond to the changes in the 
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capillary pressure as a meniscus invades into or retreats from a tube. The 
fluctuations are small compared to the total pressure and thus the capillary 
pressure cannot affect the result very much. 

In addition to containing important physics, the pressure function in 
figure H is used to establish the stability and convergence properties of the 
numerical solution. For all simulations it is verified that this function con- 
verges towards a unique solution when Ax max < O.ld and the menisci are 
updated according to the second order Runge-Kutta method as described in 
section f|. That means Kirchhoff's equation must be solved approximately 
10-20 times to let a meniscus pass through a single tube. This is proba- 
bly what we can expect when we want to measure the fluctuations in the 
pressure due to local capillary changes inside the tubes. With larger step 
lengths the variations in the capillary pressure are lost and the solution is 
no longer suitable for our measurements. 



6.2 Stable Displacement 

Figure |l0| shows the result of a simulation performed in the regime of stable 
displacement on a lattice of 60 x 60 nodes. Figure [ll| shows the corresponding 
pressure across the lattice as a function of time. The simulation is stopped 
when the invading fluid has about half filled the system. As in the case 
of viscous fingering Q = 1.5 ml/min and C a = 4.6 • 10 -3 . The viscosities 
are [i\ = 0.10 P and /U2 = 10 P giving M = 1.0 • 10 2 . Thus, the invading 
non-wetting fluid is more viscous than the defending wetting fluid. 

The fluid movements are dominated by the viscous forces in the invading 
liquid and like viscous fingering the capillary effects are negligible. The 
invading fluid generates a compact pattern with an almost flat front between 
the non-wetting and wetting fluid. The simulation is stopped after the width 
of the front has stabilized, that means after the width of the front stops 
growing. 

The average pressure across the lattice (figure |ll|) increases according 
to the amount of the high viscosity invading fluid injected into the system. 
Due to the low viscosity defending fluid the pressure corresponds to the 
pressure across the invading phase and a linear increase in the pressure is 



observed after the front has stabilized. (In figure 11: t > 50 s). Like in 
viscous fingering, the average slope of the pressure function also depends on 
the viscosity ratio M and the injection rate. 

The viscous forces dominate the pressure evolution, but fluctuations due 
to capillary effects are observed. The perturbations have about the same size 
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as for viscous fingering, which is not surprising since the size distribution of 
the radii of the tubes is the same for the two cases. The threshold pressures 
setting the strength of the capillary fluctuations, is inversely proportional 
to the radius of the tubes. 

In fast displacement the viscous forces are strong enough to deform and 
even move small regions of trapped defending fluid which are left behind 
the front. It is also observed that the invading fluid grows along the whole 
front, causing the corresponding menisci to reach the end of the tube ap- 
proximately at the same time. All together the result is a lot of bubbles of 
non-wetting and wetting fluid which travel through the network. Physically, 
the bubbles are expected but at the time they become too dominant in the 
displacement they cause some computational complication. The bubbles 
increase the number of tubes that should contain three menisci. Thus, a 
large number of reorganizations inside the tubes are applied to fulfill the 
constraint of at most two menisci in each tube. Typically, this is observed 



in figure 11 where the fluctuations in last part of the pressure function seem 
to increase in amplitude. The increasing amplitude is a result of the un- 
physical pressure jumps when three menisci are reorganized into one. For 
that reason the simulation is stopped before breakthrough of the invading 
fluid. 

Another shortcoming in fast displacement is that the effective time step 
may approach zero. When the bubbles develop the number of menisci in- 
creases causing many of them to arrive at the end of the tubes approximately 
at the same time. The result is that the difference of the arrival times for the 
menisci goes to zero and a lot of iterations are required to move the whole 
front further on. Numerical experiments have shown that the problem first 
arises when the displacement is rather fast, that means C a > 5 • 10 



-3 



6.3 The Regime of Capillary Fingering 

Finally a calibration simulation is run in the regime of capillary fingering. 



The resulting pattern is shown in figure |12j, and in figure |13| the corre- 
sponding pressure across the lattice as a function of time is plotted. The 
simulation is performed on a lattice of 40 x 60 nodes and the run is stopped 
when the invading fluid reaches the outlet. The two fluids have equal viscosi- 
ties, /ii = fj,2 = 0.50 P, corresponding to similar experiments performed by 
Frette et al. (1997). The invading non- wetting fluid displaces the defending 
wetting fluid with a rather low injection rate, Q = 0.20 ml/min. The cap- 
illary number and the viscosity ratio become C a = 4.6 • 10~ 5 and M = 1.0 
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respectively. 

In capillary fingering the displacement is so slow that the viscous forces 
are negligible, with the consequence that the main force is the capillary one 
between the two fluids. Only the strength of the threshold pressure in the 
tubes decides if the invading fluid invades that tube or not. Since the radii 
of the tubes (which determine the threshold pressures) are randomly chosen 
from a given interval, the non-wetting fluid flows a random path of least 
resistance. 

Figure 12 shows a typical rough front between the invading and the 
defending fluids with trapped clusters of defending fluid left behind the 
front. As opposed to stable displacement the clusters appear at all sizes 
between the tube length and the maximum width of the front. 

The pressure across the lattice (figure 13) exhibits sudden jumps accord- 
ing to the capillary variation when the non- wetting fluid invades (or retreats) 
a tube. The fluctuations identify the bursts where the invading fluid pro- 
ceeds abruptly. An enlargement of a small part of the pressure function at 
at time around 550 s is given in figure 0. The fi gure shows clearly this kind 
of dynamics. The pressure across the lattice slowly increases in stable peri- 
ods before the threshold pressure in the tube which is going to be invaded 
is reached. At the threshold pressure the meniscus becomes unstable and 
the invasion of fluid takes place in a burst accompanied by sudden negative 
jumps in the pressure. The pressure curve in figure [l4| is in good qualitative 
agreement with experimental results [17|. 

The main problem in capillary fingering is the computation time. The 
above simulation used about 36 CPU hours at a Cray T90, on a smaller 
lattice than the ones for viscous fingering and stable displacements. The 
two latter cases required both about 7 CPU hours. The reason for this dra- 
matic increase in computation time is explained by looking at the physical 
properties of capillary fingering. The invasion occurs in bursts and each 
burst is localized to some very few tubes. That means that even with step 
lengths equal to those of stable displacement the amount of invading fluid 
injected into the system each time step is much less. Now, compared to 
viscous fingering the total saturation of the invading fluid is relatively high 
and as a consequence, an enormous number of time steps are required to 
reach break through. In the above simulation approximately 300, 000 time 
steps are applied, while in the case of viscous fingering only about 2, 500 
steps were necessary to reach breakthrough. 
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7 Conclusions 



We have presented and discussed a two-dimensional network model simulat- 
ing drainage displacement. Moreover, we have performed numerical calibra- 
tion simulations whose results are found to be in good qualitative agreement 
with the properties observed in the three major flow regimes: viscous fin- 
gering, stable displacement and capillary fingering. 

An important feature of the model is the capability to study temporal 
evolutions of different physical properties taking place in the displacement 
process. This continues the work of understanding two-phase flow in porous 
media. In particular, we have presented calculations of the pressure evolu- 
tions where the dynamics of the capillary effects were taken into account. 
With the proposed model it is also possible to measure the viscous and 
capillary contribution to the total pressure. This is a great advance since 
corresponding measurements in experimental setups are too difficult. A 
quantification of the competition between viscous and capillary forces will 
help to characterize the different flow regimes observed in drainage displace- 
ment. 

An approximation has been developed to simulate a simultaneous flow of 
the two liquids inside a single tube. The presented moving rules take care of 
this kind of mixing and actually, they are essential to make the model work. 
So far the moving rules are based on the mechanisms observed in drainage, 
but in principle it should be possible by a certain modifications to model 
imbibition as well. 

The lattice sizes are limited by the computation time and the model 
requires usage of high performance vector machines. To afford simulations 
on lattice sizes comparable to those in experimental setups p|,|I^|, more so- 
phisticated and efficient algorithms have to be developed. Especially, vector 
machines with parallel capabilities should be considered. 
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Figure 1: A square lattice of tubes connected together at nodes. The size 
of the lattice is 10 x 10 nodes. The black region indicates the invading non- 
wetting fluid coming from below and the light gray indicates the defending 
wetting fluid flowing out of the top. 



Figure 2: Flow in a tube containing a meniscus. 



Figure 3: The capillary pressure p c as a function of the meniscus' position 
x in that tube. In the middle of the tube at x = d/2 the capillary pressure 
becomes equal to the threshold pressure pt ■ 



Figure 4: Four different fluid arrangements inside one tube. The shaded and 
the white regions indicate the non-wetting and wetting fluid respectively. 



Figure 5: The motion of the menisci at the nodes, (a): The non- wetting fluid 
(shaded) reaches the end of the tube (position 1) and is moved a distance S 
into the neighbor tubes (position 2). (b): The wetting fluid (white) reaches 
the end of the tubes (position 1) and the non-wetting fluid (shaded) retreat 
to position 2. For both (a) and (b) a proper time is recorded due to the 
small movement 5. 
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Figure 6: Reduction of three menisci into one. (a): The non- wetting fluid 
has reached the end of the tube (position 1) and is going to be moved into 
the neighbor tubes (position 2). The three menisci in the left tube are 
reorganized causing the situation shown in (b) to appear. In the figures the 
arrows denote the length of the wetting fluid between position 2 and 3 which 
is equal to the distance the meniscus at position 4 is moved to the right. 



Figure 7: A "mixture" of non-wetting (shaded) and wetting (white) which 
flow into the neighbor tubes. The different arrangements (a)-(g) are a result 
of applying the rules which are described earlier in this section. For all figures 
the fluids flow towards the node from the bottom and right tube while the 
fluids in the top and the left tube flow away from it (denoted by the arrows 
in (a)). 



Figure 8: The pattern obtained of a simulation in the regime of viscous 
fingering on a lattice of 60 x 80 nodes. C a = 4.6 • 10~ 3 and M = 1.0 ■ 10~ 3 . 
The invading non-wetting fluid (black) displaces the defending wetting fluid 
(gray) from below. The simulation took about 1\ hours on a Cray T90 
vector machine. 



Figure 9: The calculated pressure across the lattice as a function of time for 
viscous fingering. The pressure decreases with time due to the viscous forces 
in the defending fluid and the flow velocities of the moving finger tips. The 
perturbations correspond to the capillary forces due to the moving menisci. 
The time is the total time lapse required to let the invading fluid reach the 
outlet. 



21 



Figure 10: The pattern obtained of a simulation performed on a lattice of 
60 x 60 nodes in the regime of stable displacement. C a = 4.6 • 10 -3 and 
M = 1.0 ■ 10 2 . The invading non-wetting fluid (black) displacement the 
defending wetting fluid (gray) from below. The simulation took about 7 
hours on a Cray T90 vector machine. 



Figure 11: The calculated pressure across the lattice as a function of time for 
stable displacement. For t > 50 s the pressure increases linearly with time 
due to the viscous forces in the invading fluid and the constant injection 
rate. The perturbations correspond to the capillary forces along the front. 



Figure 12: The pattern obtained of a simulation performed on a lattice 
of 40 x 60 nodes in the regime of capillary fingering. C a = 4.6 ■ 10 -5 and 
M = 1.0. The invading non- wetting fluid (black) displacement the defending 
wetting fluid (gray) from below. The simulation took about 36 hours on a 
Cray T90 vector machine. 



Figure 13: The calculated pressure across the lattice as a function of time for 
capillary fingering. The fluctuations due to the capillary forces correspond 
to the burst dynamics of the invading fluid. See also figure 14. 



Figure 14: A magnification of the pressure across the lattice as a function 
of time in the interval 540-560 s. The fluid invasion takes place in bursts 
characterized by the the negative pressure jumps. 
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